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Abstract 

A parallel version of the invaded cluster algorithm is described. Results from 
large scale (up to 4096 2 and 512 3 ) simulations of the Ising model are reported. 
No evidence of critical slowing down is found for the three-dimensional Ising 
model. The magnetic exponent is estimated to be 2.482 ± .001(/3/V = 0.518 ± 
.001) for the three-dimensional Ising model. 

I. INTRODUCTION 

Advances in algorithms and hardware have greatly increased the power and scope of 
Monte Carlo studies of critical phenomena. Cluster algorithms, first introduced by Swendsen 
and Wang IJ , largely overcome the problems of critical slowing and thereby permit the study 
of very large systems using a reasonable number of Monte Carlo steps. Parallel machines 
provide platforms on which very large systems may be efficiently studied. Several groups 
0-^] have implemented cluster algorithms on parallel machines to obtain high precision 
results for the three-dimensional Ising model. 

The invaded cluster (IC) algorithm |5||J is a new type of cluster algorithm with features 
that make it very attractive for high precision studies of critical phenomena. First, it has 
less critical slowing than other cluster algorithms. Indeed, here we provide strong evidence 
that for the three-dimensional Ising model there is no critical slowing down. Second, the 
IC algorithm does not use an estimate of the critical temperature as an input; instead, the 
critical point is sampled without any fine-tuning of parameters and an estimate of the critical 
temperature is obtained as an output of the simulation. On the other hand, the IC algorithm 
samples the IC ensemble rather than the canonical ensemble. Although the two ensembles 
are believed to be equivalent in the infinite volume limit, the finite-size scaling properties 
of the IC ensemble are not yet well-understood. Here we present results that test a simple 
finite-size scaling hypothesis and determine the associated finite-size scaling exponent. 
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In this paper we describe an implementation of the invaded cluster algorithm for dis- 
tributed memory SPMD (single program multiple data) parallel machines. Using this al- 
gorithm we studied two- and three-dimensional Ising systems up to size 4096 2 and 512 3 , 
respectively. There are several motivation for this project. The first is to develop an ef- 
ficient implementation of the IC method for parallel computation. This is a non-trivial 
problem because the IC method appears to be quite sequential. The second objective is to 
better understand the dynamics of the IC algorithm and the finite-size scaling properties 
of the IC ensemble. These properties are interesting in their own right and they must be 
understood before the algorithm can be used for high-precision studies. Finally, we wish to 
obtain accurate results for the critical temperature and exponents for the three-dimensional 
Ising model. With relatively modest effort we have been able to obtain accuracy compara- 
ble to recent studies in the measurement of /3/u and somewhat less accuracy for the critical 
temperature. We argue that the IC method is likely to be the most efficient method for ob- 
taining very high precision estimates for the critical temperature and the magnetic exponent 
for the Ising model and related spin systems. 

In Sec. [TI| we describe the invaded cluster method and discuss our parallel implemen- 
tation. In Sec. [Ilj we discuss the performance of the parallel IC algorithm running on 
CM-5 machines. In Sec. [TV] we present some high-precision numerical results for the three- 
dimensional Ising model using the parallel IC algorithm. Conclusions are presented in Sec. 

m 



II. INVADED CLUSTER ALGORITHM 
A. Cluster Methods 

In cluster Monte Carlo methods a spin system is sampled in both its spin representation 
and in an associated graphical representation. The graphical representation is an interacting 
percolation model that is usually defined on the bonds of the lattice. The graphical repre- 
sentation by itself contains all of the information relevant for the study of the spin system. 
A single Monte Carlo (MC) step consists of two parts, a 'bond move' and a 'spin move.' 
During the bond move, the spins are frozen and the graphical configuration (bond config- 
uration) is updated in accord with the current spin configuration. During the spin move, 
it is the opposite: the bonds are frozen and the spin configuration is updated. In a critical 
region, tremendous improvements in the efficiency of cluster methods over traditional local 
dynamics can arise if the graphical representation has a coinciding percolation threshold. In 
this case, clusters of spins on all scales are coherently updated. 

Since the IC and Swendsen-Wang (SW) algorithms are very similar let us begin with a 
brief description of the SW algorithm. Consider the nearest neighbor Ising magnet on the 
square or cubic lattice defined by the Hamiltonian 

PH=-KJ2wj (1) 

where <7j = ±1, K is the dimensionless coupling, and the sum is over all the bonds of 
the lattice. During the bond move, each 'satisfied' bond ((i, j) is satisfied if <7j = <jj) is 
independently occupied with probability p or left vacant with probability (1 — p). Unsatisfied 
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bonds are never occupied. After a bond configuration has been created in this fashion, 
clusters are identified: A cluster is a set of (like) spins connected by occupied bonds; an 
isolated spin is also counted as clusters. Next, the spin configuration is erased and the bond 
move is completed. The spin move for the SW algorithm consists of assigning either spin-up 
or spin-down with probability 1/2 independently to every cluster. Each spin in the cluster 
takes the spin value assigned to the cluster and then the bond configuration is erased and 
the spin move is completed. The SW algorithm samples the canonical ensemble at coupling 
K where K is related to p via 

p = p(K) = 1 - e~ 2K . (2) 

The spin move for the IC algorithm is the same as for the SW algorithm: each cluster is 
assigned either spin-up or spin-down with probability 1/2. The IC bond move consists of 
occupying the satisfied bonds in random order until the cluster configuration first fulfills a 
stopping rule. During the bond move, bonds are occupied one at a time and after each bond 
is added, the cluster structure is checked to see if the stopping rule is satisfied. There are 
a number of useful stopping rules, some of which, evidently, cause the algorithm to sample 
the Ising critical point. In this study, we will use 'topological' rules which are appropriate 
for systems with periodic boundary conditions. There are several versions: The '1-span' 
topological rule demands that some cluster is connected around the system in at least one 
of the periodic directions - as soon as this condition is fulfilled, the bond move is complete 
and no more bonds are added to the configuration. In general, the fc-span rule (k < d) 
requires that there is connectivity around the system in at least k directions. Although all 
of the fc-span rules sample the Ising critical point, the different rules have different finite-size 
scaling behavior. 

The critical coupling and the magnetic exponent may be obtained in the IC approach. 
We begin with K c . Let / be the ratio of the number bonds occupied during the bond move 
to the number of satisfied bonds for a single Monte Carlo (MC) step. The average value of 
/ is an estimator of the critical temperature according to the relation 

lim </>= 1 -e~ 2Kc (3) 

L^oo 

where K c is the critical coupling and L is the system size. To understand this relation, 
note that the only difference between the SW and IC algorithms is the way in which / is 
determined in a bond move. For the SW algorithm, each satisfied bond is independently 
occupied with probability p so that the ratio of occupied to satisfied bonds is a binomial 
random variable. For the IC algorithm, the stopping rule determines the number of occupied 
bonds. However, for large systems, one can presume that the fluctuations in / become small 
for both algorithms. As a consequence, the IC algorithm that outputs < f> yields the same 
value for all local observables as the SW algorithm with p{K) =< f > so that the relation 
< / >~ 1 — e~ 2K yields an estimate of the temperature of the IC ensemble. The bond 
percolation threshold on the set of satisfied bonds at the Ising critical point is known to be 
p(K c ). Since the /c-span rules forces the cluster configuration to the percolation threshold, 
we see that Eq. @ holds. This argument (which still lacks a rigorous proof) is discussed in 
more detail in Ref. || 

Let us now turn to a discussion of the magnetic exponent. For a critical stopping rule, 
such as a topological rule, the cluster that fulfills the stopping condition is called the spanning 
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cluster. This spanning cluster can be thought of as a typical example of a large-scale cluster 
at criticality and thus its fractal dimension yields the magnetic exponent for the system. 
In particular, if M is the number of sites in the spanning cluster, we determine the fractal 
dimension D via 

M ~ L D (4) 

The fractal dimension is related to other critical exponents via D = y^ = d — (3/v where y^ 
is the magnetic exponent, (3 the magnetization exponent, v the correlation length exponent 
and d the spatial dimension. At the present time we do not know how to extract the thermal 
exponent from the IC algorithm. 

The finite-size scaling properties of the IC ensemble are currently under investigation. 
Here we state the fundamental assumption that will be used later in the analysis of the data. 
Let J-(f, L) be the (cumulative) distribution function for / measured in a system of size L. 
We suppose that there is an exponent, u and function Q such that 

F{f,L)~Q[{f-p{K c ))L u ] (5) 

Naive finite-size scaling arguments suggest that u= 1/v with v the Ising correlation length 
exponent, however it is quite apparent that this is not the case. Measurements of u reported 
in Ref. B and below in Sec. ffVI yield a value u ~ 0.69 in d = 3. Presumably this u is the 



inverse of the correlation length exponent for a related percolation-Ising hybrid model but 
at present, no definitive conclusions have been obtained. 



B. Parallelizing the Invaded Cluster Algorithm 

The approach taken here to parallelize the IC algorithm is based upon the parallel SW 
algorithm described in Ref. ||. The basic idea for both algorithms is to divide the lattice 
into cells each of which is handled by a single processor. During the bond move the cluster 
configuration must be determined. Every site of the lattice has an integer label that identifies 
its cluster membership. Two sites are in the same cluster if they have the same label. At 
the beginning of the bond move, every site has a unique local label. The SW bond move 
proceeds in two stages. In the local stage, each processor in parallel occupies a fraction p of 
the satisfied bonds in its cell and identifies the local cluster structure using a conventional 
sequential method. The local clusters are the clusters that would exist if no bonds on the 
boundaries of the cell were occupied. During the relaxation stage, global labels are created 
and global clusters are identified. Global clusters which are formed out of local clusters 
by boundary connections between cells are assigned the lowest valued label from the the 
labels of the connected local clusters. This is accomplished iteratively by the exchange of 
labels between neighboring cells. A very similar relaxation stage is used for the IC algorithm 
except that more information must be exchanged between nodes to allow for checking the 
stopping condition. 

Since communication between processing nodes is much slower than processes handled 
by a single node it is desirable that individual nodes do as much work as possible between 
communication steps. On the other hand, the sequential IC algorithm checks for spanning 
after each new bond is occupied. Clearly this would not be a good approach to directly 
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parallelize. Instead, the algorithm converges to the value of / by a sequence of choices of 
p. Imagine that an independent uniform random number between and 1 is assigned to 
each bond of the lattice. (In fact, this is not quite the approach that is actually used but, 
for purposes of explanation, suppose that this is the case.) On the first step each processor 
provisionally occupies all satisfied bonds with random numbers less than p\. The resulting 
cluster configuration is identified and the stopping rule is checked. For the sake of argument, 
suppose that spanning is not detected. In this case, the provisionally occupied bonds are 
permanently occupied and a new value p<i > p\ is chosen. The current cluster configuration 
is stored and additional bonds are provisionally occupied with random numbers from pi to 
Pi- The cluster configuration is modified accordingly and the stopping condition is checked. 
If spanning is not detected, the new cluster configuration is stored and the provisionally 
occupied bonds given tenure. The value of p is incremented in this way until on the sth step 
the stopping condition is fulfilled. At this point we know that spanning occurs between p s -\ 
and p s . We backtrack to the cluster configuration (saved in memory) associated with p s -\ 
and try again, this time choosing p = (p s -i +p s )/2. Henceforth, the algorithm carries out a 
binary search, each time reducing the range of p by half until the stopping condition is just 
fulfilled and / is determined. 

In the above discussion we supposed that random numbers are assigned to each bond 
and that on each step all satisfied bonds with random numbers less than p are occupied. 
This is not an efficient procedure because it requires examining every bond each time p is 
changed. Instead we do the following. Each processor creates a random permutation of the 
n bonds in its cell. This permutation determines the order in which bonds will be tested 
by the processor. On the first step, processor i tests its first Kp bonds where is chosen 
from the binomial distribution, B[n,pi] for n trials with success probability p\. This is 
statistically equivalent to testing all bonds with random numbers less than p 1 as described 
in the previous paragraph. Suppose that spanning is not found, then p\ is incremented to p2- 
At this point processor i has permanently occupied op ^— h[ bonds. Processor i now tests 
the next hf* bonds where now h[ is chosen from B[n — ,{p2—pi)/ (1 — Pi)] • This process 
of occupying bonds continues as long as the stopping condition is not fulfilled and on the mth 
such step the number of occupied bonds is incremented according to <— a- m ~ 1 ' ) + h\ m \ 
Suppose that spanning is first detected on step s and that h\ bonds have been provisionally 
occupied by processor i during this step. On the next step each processor tries to occupy a 
smaller number of bonds corresponding to splitting the difference between p s -i and p s . To 
do this, processor i tests the next h!f +l ^ bonds starting with the (a\ s ^ + l)th bond where 
h\ s is chosen from B\h\ , 1/2]. A little reflection shows that this procedure is statistically 
equivalent to the one described above for carrying out the binary search. The advantage of 
this approach is that on successive steps it is not necessary to check every bond in the cell. 
Note that the local ordering of bonds is fully specified at the outset of the bond move but 
that the global ordering of the bonds is generated as needed and never fully specified. 

III. PERFORMANCE 

We tested the parallel IC algorithm for the two-dimensional Ising model on a 32-node 
CM-5 machine. Fig. [I] shows that the number of relaxation cycles n re ia X per MC step as a 
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function of system size L. The roughly logarithmic increase in n re i ax results from the binary 
search method used to find /. We separately measured the time spent in local cluster growth 
tiocai and relaxation t re i ax - The local time increases as the number of spins per node, thus, for a 
fixed number of processors, £i oca i ~ L d . The relaxation time is dominated by communication 
between nodes so t re iax is proportional to the boundary area of the local lattice and the 
number of relaxation cycles, t ro i a x ~ L d ~ x logL. Fig. § shows local and relaxation times as 
a function of system size. The relaxation time is greater for small systems but, when the 
number of spins per node exceeds several hundred thousand the local time dominates. 

In Fig. |3] we compare the performance of the parallel IC algorithm with that of the 
parallel SW algorithm for the two-dimensional Ising model. The total time per spin per 
MC step for both algorithms decreases as the system size increases and saturates for large 
systems. For large systems in both two and three dimensions, the time for one MC step for 
the parallel IC algorithm is about a factor of 7 larger than for the parallel SW algorithm. 

We obtained 1.1 /i s per spin per MC step for the parallel IC algorithm for the two- 
dimensional Ising models with system size 4096 2 on the 32-node CM-5 machine and 3.4 \x s 
per spin per MC step for the three-dimensional Ising model on the same machine and with 
the same number of spins (256 3 ). The slower speed for three-dimensions results from the 
need for more communications between nodes since each node has six neighbors rather than 
four. For both two and three dimensions the time per spin per MC step decreases as the 
system size increases. The relaxation time is, however, still comparable to the local time for 
256 3 , so the parallel IC algorithm for the three-dimensional Ising model can be expected to 
become more efficient as the system size increases holding the number of processors fixed. 

To test the sensitivity of the algorithm to the random number generator, we performed 
simulations using two different random number generators. Most of the data is collected 
using the 250-shift generator in which new random numbers are produced by Xi = Aj_ 10 3 © 
-Xj-250) where Xi are 32-bit integers and © represents the bitwise exclusive OR operation. 
We also tested a combined shift register generator, introduced in Ref. ||. The combined 
generator consists of two shift registers dj = cii-9218 © and 6, = 61-97 © h-m with 

32-bit integers a*, b^. New random numbers of the generator are combined by Xi = 0,4 © fcj. 
We observed no difference between results from the two generators. For example, the mean 
value of / for system size 80 3 was measured as .358 144 ± .000 006 for the first generator 
and .358 143 ± .000 003 for the second. 

IV. NUMERICAL RESULTS 

We studied the three-dimensional Ising model using the parallel IC algorithm for system 
sizes 16 3 to 512 3 . Unless otherwise stated we used the 1-span rule. For each run we allowed 
200 MC steps for equilibration before collecting data. Independent runs consisted of 10 000 
MC steps for the smallest systems (up to 96 3 ) to 500 MC steps for the largest system (512 3 ). 
For each size, a number of independent runs were made. We obtained statistics for / and 
the size of the spanning cluster. Error bars represent the standard error (one standard 
deviation) associated with these independent runs. The data for the 1-span rule is given in 
Table I. 

Figure (| is a double logarithmic plot of aj = (< f 2 > — < f > 2 ) 1 ! 2 vs. system size. A 
fit of the form a f = a + bL- u , yields a = 0.0003 ± .0002 and u = 0.69 ± .01. This result 
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strongly supports the hypothesis that Uf approaches zero as a power law. 

The finite-size scaling hypothesis, Eq. (|5]) implies an asymptotically linear relation be- 
tween < f > and (j j. Figure || shows plots of < / > against <T/ for the three /c-span rules. 
Although these plots are consistent with asymptotic linear behavior, it is clear that there 
are large subleading corrections to scaling. For the 1-span rule these corrections are negative 
and for the 3-span rule they are positive. 

Figure ^| shows more detail for the 1-span rule where, over the range of system sizes 
studied, < f> varies only in the fifth significant digit. However, to obtain yet higher accuracy 
than p(K c ) = .3581 requires either larger system sizes or a good fit to the subleading terms. 
Fitting to the form <f>= p(K c ) + aa f + baj we obtain p(K c ) = 0.358 068 ± .000 009 with 
a = .021 and b = —1.47. The corresponding critical coupling, quoted with two standard 
deviation errors, is K c = 0.221 637 ± .000016. The fit is shown as a dotted line in the 
figure. The primary difficulty in obtaining a high precision value of the critical coupling is 
the relatively large corrections to scaling, represented by the fact that the fitting parameter 
b is a good deal larger than a. Our estimate of the critical coupling is consistent with recent 
values (e.g. K c = 0.221 654 6 ± .000 001 §) but the error is much larger. 

Figure [7] shows the median value of / for the 1-span rule versus o~/. According to the 
finite-size scaling hypothesis this plot should also be linear and, indeed, approximately linear 
behavior is evidenced over the whole range of system sizes however it should be noted that 
the variation in the median is much larger than for the mean. Furthermore, errors in the 
mean are less. Thus appears that the mean is a better estimator of the critical coupling. 

We measured the mass (the number of sites) M of the spanning clusters. Using Eq. (|J) 
we can find the fractal dimension of the spanning cluster and the magnetic exponent. The 
double logarithmic plot of M and L in Fig. ||| shows that the data is well described by a 
power law and an examination of different ranges of system size yields no clear trends in 
the power. A linear fit for L > 48 yields (3/v = 0.518 ± .001 consistent with a recent value 
0.5185 ± .0008 §. 

Next we consider the dynamical properties of the algorithm. The normalized autocorre- 
lation function of an observable A is defined by, 

<A A t > - <A> 2 
TA{t) = <A*>-<A>* ' (6) 

where t is time in Monte Carlo steps. Figure |9| shows the autocorrelation function for / 
for a 128 3 system. The striking feature of this curve is that it is negative on time step 
1 indicating anti-correlations between successive MC steps. As discussed in Refs. P,pL0| 
these anti-correlations are associated with the negative feedback mechanism that forces IC 
dynamics to the critical point. By the fourth step it is not possible to distinguish the 
autocorrelation function from zero. 

The integrated autocorrelation time ta is defined by, 

1 w 

r A (w) = - + ^A(t), (7a) 
A t=i 

r A = lim t a (w). (7b) 
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Using the procedure described in Ref . , we measured integrated autocorrelation times 
for the magnetization m, energy e and temperature estimator /. The autocorrelation times 
are given in Table II and plotted against system size in Fig. [10]. Due to these anti-correlations 
Tf and t 6 decrease with system size; r m remains nearly constant. Thus IC dynamics for the 
three-dimensional Ising model apparently does not suffer critical slowing. (Although, as 
is the case for all cluster algorithms, the parallel time needed for a single MC step must 
increases logarithmically with the system size [ID]). The dynamic exponent z is defined 



by r ~ L z . We obtained the dynamic exponents z m = —.01 ± .01, z e = — .46 ± .01 and 
Zf = — 1.23 ± .10 for m, e and / respectively. 

Figure [11] shows Tj(l) as a function of 1/L. It appears that the infinite volume limit 
for Tf(l) lies between between —0.6 and —0.65. Since Tf > there must be corresponding 
positive values summing to at least 0.1 to 0.15 for t > 2. It would be interesting to have 
more information about the infinite volume limit of the / autocorrelation function. 

The autocorrelation function is related to the standard error, 8 A in measuring < A > 
according to, 

-r oo 

8 A = a A {2/Nfl\r A - ^£*r A (i)) 1/2 (8) 

with a\ =< A 2 > — <A> 2 and iV the number of Monte Carlo steps. In most Monte Carlo 
applications the second term is much smaller than the first however for the IC algorithm this 
may not be the case since r/, r e — > 0. The autocorrelation function for / is approximately 
— 1/2 for t — 1 and small for t > 1 so that the approximate error in measuring the critical 
coupling is 

8f « a f {2/Nf'\r f + ±)W (9) 

Suppose that the number of MC steps is large enough that Tf ^> 1/N as is the case here 
for the smaller system sizes. Combining the scaling behavior of Tf and of Gf we have, for 

8f ~ L {z i /2 - u) /N 1/2 ~ L- L31 /iV 1/2 . (10) 

On the other hand, if Tf and 1/N are comparable as is the case here for the largest lattices 
then 5f ~ L^ z f- U ^ = L~ im and N ~ L~ z f = L L23 . 



V. CONCLUSIONS 

We used a parallel implementation of the invaded cluster algorithm to study two- and 
three-dimensional Ising models. We have found that the parallel algorithm is quite efficient 
for large system sizes where most computational work is carried out constructing clusters 
within processing nodes and relatively little time is spent in communication. 

We have shown that the invaded cluster algorithm has no critical slowing down for the 
three-dimensional Ising model as measured by integrated autocorrelation times for the mag- 
netization, energy or estimated critical temperature. Indeed, the integrated autocorrelation 
time for the latter two quantities approaches zero as the system size grows so that relatively 
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few Monte Carlo steps are needed for high accuracy estimates of the finite-size critical energy 
or critical temperature. 

On the other hand, extrapolation to the infinite volume limit for quantities such as the 
critical temperature is complicated by the absence of a well-understood finite-size scaling 
theory for the invaded cluster ensemble. We have examined the hypothesis that the distri- 
bution of / is controlled by a single scale factor and found that this apparently holds for 
large systems (greater than 128 3 ) but that there are significant corrections to the predicted 
asymptotic linear dependence of <f> on aj. This is not unexpected since the width of the 
/ distribution is rather large compared to | <f> —p(K c )\. Thus relatively small corrections 
to scaling for the / distribution can lead to relatively large deviations from the proposed 
linear behavior of <f> vs. <7/. This means that to confidently reduce systematic errors in 
estimating the infinite volume critical temperature we would either have to go to much larger 
systems to make Of small or we would need a detailed understanding of the corrections to 
scaling in the IC ensemble. Our current estimates of the critical temperature using the IC 
method are consistent with the best conventional estimates but have errors which are about 
an order of magnitude larger. 

In contrast to the critical temperature, corrections to scaling in measuring the critical 
dimension of the spanning cluster are very small and we are able to obtain an excellent 
estimate of the magnetic exponent. Our result is y h = 2.482 ± .001(/3/V = 0.518 ± .001). 
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FIG. 1. The number of relaxation cycles per MC step n re i a x vs. system size L for the 
two-dimensional Ising model. 
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FIG. 3. Time per spin per MC step for the two-dimensional Ising model vs. L for the parallel 
IC and SW algorithms. 
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FIG. 4. Double logarithmic plot of 07 vs. system size L for the three-dimensional Ising model. 
The slope of the curve is .69 ± .01 from a linear fit. 
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FIG. 5. Plot of / vs. Of for three fc-span rules for the three-dimensional Ising model. The point 
on the vertical axis is the accepted infinite volume estimate. 
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FIG. 6. <f> vs. Of for the three-dimensional Ising model using the 1-span rule. The dotted 
line is a quadratic fit and yields p{K c ) = 0.358 068 ± .000009. The point on the vertical axis is the 
accepted infinite volume estimate. 
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FIG. 7. The median value of / vs. cry for the three-dimensional Ising model using the 1-span 
rule. 
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FIG. 8. Double logarithmic plot of the mass of the spanning cluster, M vs. system size L for 
the three-dimensional Ising model. 
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FIG. 9. Autocorrelation function for f,Tf vs. t for the three-dimensional Ising model for 
system size 128 3 . 
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FIG. 10. Double logarithmic plot of autocorrelation times r m ,r e ,r/ vs. system size L for the 
three-dimensional Ising model. 
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FIG. 11. Plot of Tj(l) vs. L 1 for the three-dimensional Ising model. 
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TABLES 



TABLE I. Numerical data for the three-dimensional Ising model using the 1-span rule. 



T 

L 




if) 


r 

/median 


TV T 

M 


e 


16 


.02645(4) 


.357367(15) 


.356223(33) 


859(2) 


-1.995544(57) 


32 


01 599(3) 


357961 (11) 


357205(33) 


4830(4) 


-1 995055(26) 


48 


.01206(2) 


.358106(5) 


.357534(15) 


13256(9) 


-1.995242(23) 


64 


.00988(2) 


.358128(4) 


.357621(14) 


27078(25) 


-1.995289(25) 


80 


.00855(2) 


.358145(3) 


.357696(11) 


47019(43) 


-1.995339(19) 


96 


.00737(2) 


.358144(5) 


.357783(17) 


74036(86) 


-1.995362(33) 


128 


.00625(3) 


.358140(2) 


.357851(16) 


151241(208) 


-1.995377(18) 


160 


.00540(2) 


.358136(3) 


.357898(19) 


263220(505) 


-i.yyoooy^io j 


256 


.00400(4) 


.358130(2) 


.357904(21) 


847210(3238) 


-1.995425(25) 


512 


.00250(7) 


.358105(11) 


.358098(36) 


4668684(66195) 


-1.995340(40) 


TABLE II. Integrated autocorrelation 


times for three-dimensional Ising 


models for the IC 


algorithm. Results 


are measured at time step w = 6. 






L 












16 




.325(4) 




.639(3) 


.097(4) 


32 




.248(5) 




.637(5) 


.048(6) 


48 




.201(4) 




.628(4) 


.028(5) 


64 




.173(5) 




.636(5) 


.018(5) 


80 




.153(5) 




.630(5) 


.016(6) 


128 




.107(8) 




.618(8) 


.008(8) 


160 




.094(11) 




.606(11) 


.008(12) 
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